Square Wave Oscillator
System type: linear hybrid
State dimension: 1
Application domain: Electronics
Model description
Hybrid automaton formulation
using ReachabilityAnalysis, Symbolics, Plots
LazySets.set_ztol(Float64, 1e-15)
function multistable_oscillator(; X0 = Interval(0.0, 0.05),
V₊ = +13.5, V₋ = -13.5,
R = 20.E3, C = 5.5556E-8,
R1 = 20.E3, R2 = 20.E3)
@variables x
τ = 1/(R*C)
α = R2/(R1+R2)
A = -τ
b = (τ/α) * V₊
I₊ = HalfSpace(x <= α*V₊)
m1 = @system(x' = Ax + b, x ∈ I₊)
b = (τ/α) * V₋
I₋ = HalfSpace(x >= α*V₋)
m2 = @system(x' = Ax + b, x ∈ I₋)
automaton = LightAutomaton(2)
add_transition!(automaton, 1, 2, 1)
g1 = Hyperplane(x == α*V₊)
r1 = ConstrainedIdentityMap(1, g1)
add_transition!(automaton, 2, 1, 2)
g2 = Hyperplane(x == α*V₋)
r2 = ConstrainedIdentityMap(1, g2)
modes = [m1, m2]
resetmaps = [r1, r2]
H = HybridSystem(automaton, modes, resetmaps)
# initial condition in mode 1
X0e = [(1, X0)]
return IVP(H, X0e)
endmultistable_oscillator (generic function with 1 method)
Results
prob = multistable_oscillator()
sol = solve(prob, T=100e-4, alg=INT(δ=1.E-6), fixpoint_check=false);
plot(sol, vars=(0, 1), xlab="t", ylab="v-")
tspan.(sol)19-element Array{IntervalArithmetic.Interval{Float64},1}:
[0, 0.000320001]
[0.000316999, 0.000888001]
[0.000883999, 0.00145601]
[0.00145099, 0.00202401]
[0.00201799, 0.00259201]
[0.00258499, 0.00316001]
[0.00315199, 0.00372801]
[0.00371899, 0.00429601]
[0.00428599, 0.00486401]
[0.00485299, 0.00543201]
[0.00541999, 0.00600001]
[0.00598699, 0.00656801]
[0.00655399, 0.00713601]
[0.00712099, 0.00770401]
[0.00768799, 0.00827201]
[0.00825499, 0.00884001]
[0.00882199, 0.00940801]
[0.00938899, 0.00997601]
[0.00995599, 0.0100211]location.(sol)19-element Array{Int64,1}:
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1Let us analyze in some detail the first transition. If we plot the last 10 reach-sets of the first flowpipe, we observe that only the last 3 actually intersect the guard:
plot(sol[1][end-10:end], vars=(0, 1), xlab="t", ylab="v-")
plot!(x -> 6.75, xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red)
We now cluster those reach-sets into a single hyperrectangle:
Xc = cluster(sol[1], [318, 319, 320], BoxClustering(1))1-element Array{ReachSet{Float64,Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}},1}:
ReachSet{Float64,Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}}(Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}([6.7258286982437765], [0.024171301756223507]), [0.000316999, 0.000320001])Plotting Xc matches with the flowpipe after the jump:
plot(sol[1][end-10:end], vars=(0, 1))
plot!(sol[2][1:10], vars=(0, 1))
plot!(x -> 6.75, xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red)
plot!(Xc[1], vars=(0, 1), c=:grey)
Finally, we note that the algorithm finds an invariant of the system after the first period. To activate such check pass the fixpoint_check=true flag to the hybrid solve API.
sol = solve(prob, T=100e-4, alg=INT(δ=1.E-6), fixpoint_check=true);plot(sol, vars=(0, 1), xlab="t", ylab="v-")
When the fixpoint check is activated, the computation terminates as soon as the last reach-set is contained in a previously explored initial state.
tspan(sol)[0, 0.00145601]